#ifndef ND_POLYNOM_H
#define ND_POLYNOM_H

#include <vector>
#include <algorithm>
#include <stdlib.h>
#include <QDebug>

#define EPS 1e-8

template< class C, int N >
class nd_term
{
	C       m_coeff;
	int     m_pow[N];
public:
	nd_term()
	{
		m_coeff = 0;
		for( int i = 0; i != N; ++i )
			m_pow[i] = 0;
	}

	nd_term( const nd_term& other )
	{
		m_coeff = other.m_coeff;
		for( int i = 0; i != N; ++i )
			m_pow[i] = other.m_pow[i];
	}

	nd_term( const C& coeff, const int* pows )
	{
		m_coeff = coeff;
		for( int i = 0; i != N; ++i )
			m_pow[i] = pows[i];
	}

	bool operator < ( const nd_term& other ) const
	{
		for( int i = 0; i != N; ++i )
		{
			if( m_pow[i] > other.m_pow[i] )
				return true;
			else if( m_pow[i] < other.m_pow[i] )
				return false;
		}
		return m_coeff < other.m_coeff;
	}

	bool is_const() const
	{
		for( int i = 0; i != N; ++i )
		{
			if( m_pow[i] != 0 )
				return false;
		}
		return true;
	}

	bool is_zero( const C& eps ) const
	{
		return fabs( m_coeff ) <= eps;
	}

	void mul( const nd_term& other )
	{
		for( int i = 0; i != N; ++i )
		{
			m_pow[i] += other.m_pow[i];
		}
		m_coeff *= other.m_coeff;
	}

	bool is_divide_by( const nd_term& other) const
	{
		for( int i = 0; i != N; ++i )
		{
			if( m_pow[i] < other.m_pow[i] )
				return false;
		}
		return true;
	}

	bool div( const nd_term& other )
	{
		if( !is_divide_by( other ) )
			return false;
		for( int i = 0; i != N; ++i )
		{
			m_pow[i] -= other.m_pow[i];
		}
		m_coeff /= other.m_coeff;
		return true;
	}

	static bool resolve_engagement( const nd_term& m1, const nd_term& m2, nd_term& k1, nd_term& k2 )
	{
		k1.m_coeff = m2.m_coeff;
		k2.m_coeff = m1.m_coeff;

		for( int i = 0; i != N; ++i )
		{
			k1.m_pow[i] = std::max( m2.m_pow[i] - m1.m_pow[i], 0 );
			k2.m_pow[i] = std::max( m1.m_pow[i] - m2.m_pow[i], 0 );
//			qDebug() <<	"m1" <<	m1.m_pow[i] << "m2" <<	m2.m_pow[i] << "k1" <<	k1.m_pow[i] << "k2" <<	k2.m_pow[i];
		}

		for( int i = 0; i != N; ++i )
		{
			if( k1.m_pow[i] != m2.m_pow[i] ||
				k2.m_pow[i] != m1.m_pow[i] )
				return true; //Have engagement
		}
		return false; //No engagement
	}

	bool is_pow_eq( const nd_term& other ) const
	{
		for( int i = 0; i != N; ++i )
		{
			if( m_pow[i] != other.m_pow[i] )
				return false;
		}
		return true;
	}

	C get_coeff() const
	{
		return m_coeff;
	}

	int get_pow( int i ) const
	{
		return m_pow[i];
	}

	void set_coeff( const C& coeff )
	{
		m_coeff = coeff;
	}

	template< class Stream >
	void print( Stream s ) const
	{
		static const char*	names[] = { "x", "y", "z", "u", "v", "w", "t" };
		s << m_coeff;
		for( int i = 0; i != N; ++i )
		{
			int	pow = get_pow( i );
			if( pow > 0)
				s << "*";
			if( pow == 0 )
			{
			}
			else
			if( pow == 1 )
			{
				s << names[i];
			}
			else
			{
				s << names[i] << "^" << pow;
			}
		}
	}
};

template< class C, int N >
class nd_polynom
{
	std::vector< nd_term<C,N> > p;

public:

	bool is_empty() const
	{
		return p.size() == 0;
	}

	void pop_front()
	{
		if( p.size() > 0 )
			p.erase( p.begin() );
	}

	bool is_const() const
	{
		if( p.size() == 0 )
			return true;
		for( size_t n = 0; n != p.size(); ++n )
		{
			if( !p[n].is_const() )
				return false;
		}
		return true;
	}

	void clear()
	{
		p.clear();
	}

	size_t size()
	{
		return p.size();
	}

	const nd_term<C,N>& superior() const
	{
		return p.front();
	}

	void remove_superior()
	{
		p.erase( p.begin() );
	}

	void add( const nd_term<C,N>& t )
	{
		for( size_t n = 0; n != p.size(); ++n )
		{
			nd_term<C,N>&	pt = p[n];
			if( pt.is_pow_eq( t ) )
			{
				pt.set_coeff( pt.get_coeff() + t.get_coeff() );
				return;
			}
		}
		p.push_back( t );
	}

	void sub( const nd_term<C,N>& t )
	{
		for( size_t n = 0; n != p.size(); ++n )
		{
			nd_term<C,N>&	pt = p[n];
			if( pt.is_pow_eq( t ) )
			{
				pt.set_coeff( pt.get_coeff() - t.get_coeff() );
				return;
			}
		}
		p.push_back( t );
		p.back().set_coeff( -t.get_coeff() );
	}

	void mul( const nd_term<C,N>& t )
	{
		for( size_t n = 0; n != p.size(); ++n )
		{
			nd_term<C,N>&	pt = p[n];
			pt.mul( t );
		}
	}

	bool div( const nd_term<C,N>& t )
	{
		for( size_t n = 0; n != p.size(); ++n )
		{
			nd_term<C,N>&	pt = p[n];
			if( !pt.div( t ) )
				return false;
		}
		return true;
	}

	void add( const nd_polynom& s )
	{
		for( size_t n = 0; n != s.p.size(); ++n )
		{
			const nd_term<C,N>&	pt = s.p[n];
			add( pt );
		}
	}

	void sub( const nd_polynom& s )
	{
		for( size_t n = 0; n != s.p.size(); ++n )
		{
			const nd_term<C,N>&	pt = s.p[n];
			sub( pt );
		}
	}

	void chop( const C& eps )
	{
		if( p.size() == 0 )
			return;
		for( size_t n = p.size(); n > 0;  )
		{
			--n;
			if( p[n].is_zero( eps ) )
				p.erase( p.begin() + n );
		}
		if( p.size() > 1 )
		{
			std::sort( p.begin(), p.end() );
		}
	}

	static bool spolynom( const nd_polynom& f1, const nd_polynom& f2, nd_polynom& s )
	{
		if( f1.p.size() == 0 )
			return false;
		if( f2.p.size() == 0 )
			return false;

		const nd_term<C,N>&		m1 = f1.superior();
		const nd_term<C,N>&		m2 = f2.superior();
		nd_term<C,N>			k1,k2;

		if( nd_term<C,N>::resolve_engagement( m1, m2, k1, k2 ) )
		{
			nd_polynom	tmp( f2 );
			s = f1;

			//! remove superior monom
			tmp.pop_front();
			s.pop_front();

			tmp.mul( k2 );
			s.mul( k1 );

			s.sub( tmp );
			s.chop( EPS );
			return true;
		}
		return false;
	}

	template< class Stream >
	void print( Stream s ) const
	{
		static const char*	names[] = { "x", "y", "z", "u", "v", "w", "t" };
		for( size_t n = 0; n != p.size(); ++n )
		{
			bool	need_coeff = fabs( fabs( p[n].get_coeff() ) - 1 ) > 1e-16;
			bool	need_star = need_coeff;

			if( p[n].get_coeff() > 0 )
			{
				s << "+";
				if( need_coeff )
					s << fabs( p[n].get_coeff() );
			}
			else
			if( p[n].get_coeff() < 0 )
			{
				s << "-";
				if( need_coeff )
					s << fabs( p[n].get_coeff() );
			}

			for( int i = 0; i != N; ++i )
			{
				int	pow = p[n].get_pow( i );
				if( need_star && pow > 0)
					s << "*";
				if( pow == 0 )
				{
				}
				else
				if( pow == 1 )
				{
					s << names[i];
				}
				else
				{
					s << names[i] << "^" << pow;
				}
				if( pow > 0 )
					need_star = true;
			}
			if( !need_star )
				s << 1;
		}
	}
};

template< class C, int N >
class nd_polynom_system
{
	std::vector< nd_polynom<C,N> >	polynomials;
public:
	void push_back( const nd_polynom<C,N>& p )
	{
		polynomials.push_back( p );
	}

	void reduce( const nd_polynom<C,N> f, nd_polynom<C,N>& r ) const
	{
		nd_polynom<C,N>		p(f);
		size_t				i;

		r.clear();

		while( !p.is_empty() )
		{
			for( i = 0; i != polynomials.size(); ++i )
			{
				const nd_polynom<C,N>	curr = polynomials[i];

				if( curr.is_empty() )
					continue;

				const nd_term<C,N>&		f_sup = curr.superior();
				const nd_term<C,N>&		p_sup = p.superior();

				if( p_sup.is_divide_by( f_sup ) )
				{
					nd_term<C,N>		m( p_sup );
					nd_polynom<C,N>		tmp( curr );
					m.div( f_sup );

					tmp.mul( m );
					p.sub( tmp );
					p.chop( EPS );
					break;
				}
			}
			if( i == polynomials.size() )
			{
				r.add( p.superior() );
				p.remove_superior();
			}
		}
		r.chop( EPS );
	}

	void construct_groebner_basis()
	{//Buchberger
		if( polynomials.size() < 2 )
			return;

		for( size_t i = 0; i < polynomials.size(); ++i )
		{
			for( size_t j = i + 1; j < polynomials.size(); ++j )
			{
				const nd_polynom<C,N>&	pi = polynomials[i];
				const nd_polynom<C,N>&	pj = polynomials[j];

				nd_polynom<C,N>			s;
				nd_polynom<C,N>			r;

				if( !nd_polynom<C,N>::spolynom( pi, pj, s ) )
					continue;
				reduce( s, r );

				if( !r.is_empty() )
				{
					push_back( r );
				}
			}
		}
	}

	void minimize_groebner_basis()
	{
		if( polynomials.size() < 2 )
			return;

		size_t pass_count = polynomials.size();
		for( size_t pass = 0; pass < pass_count; ++pass )
		{
			bool		erased = false;
			for( size_t i = 0; i < polynomials.size(); ++i )
			{
				const nd_term<C,N>&		i_sup = polynomials[i].superior();
				erased = false;
				for( size_t j = 0; j < polynomials.size(); ++j )
				{
					if( i == j )
						continue;
					const nd_term<C,N>&		j_sup = polynomials[j].superior();
					if( j_sup.is_divide_by( i_sup ) )
					{
						polynomials.erase( polynomials.begin() + j );
						erased = true;
						break;
					}
				}
				if( erased )
					break;
			}
			if( !erased )
				break;
		}
	}

	template< class Stream >
	void print( Stream s ) const
	{
		s << "";
		for( size_t n = 0; n != polynomials.size(); ++n )
		{
			s << n << ") ";
			polynomials[n].print( s );
			s << "\r\n";
		}
	}
};

#endif // ND_POLYNOM_H
